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ABSTRACT 

We explore galaxy properties and their link with environment and clustering using a 
population of 1000 galaxies formed in a high resolution hydrodynamic simulation of 
the ACDM cosmology. At the redshift we concentrate on, z = 1, the spatial resolution 
is 1.4 proper /i~^kpc and the mass resolution is such that Milky-way sized disk galaxies 
contain ^ 10^ particles within their virial radii. The simulations include supermassive 
black hole accretion and feedback as well as a multiphase model for star formation. 
Overall, we find that a number of familiar qualitative relationships hold approximately 
between galaxy properties, for example we observe galaxies as lying between two 
broad extremes of type, where "late" types tend to be smaller in size, have lower 
circular velocities, younger stars, higher star formation rates, larger disk to bulge 
ratios and lower Sersic indices than "early types" . We find that as in previous studies 
of small numbers of resimulated galaxies the stellar component of disk galaxies is not 
as rotationally supported as in observations. Bulges contain too much of the stellar 
mass, although exponential disks do have scale lengths compatible with observations, 
and the stellar mass TuUy-Fisher relation at z = 1 is reproduced. The addition of 
black hole physics to the simulations does not appear to have made a large impact 
on the angular momentum results compared to these other studies, nor do we find 
that running an identical simulation with significantly lower mass resolution affects 
this aspect. Despite this, we can profitably use the rank order of either disk to total 
ratio, Sersic index, or galaxy age to separate galaxies into morphological classes and 
examine the density-morphology relation and morphology dependence of clustering. 
We find that while at redshift = 0, the well known preponderance of early types in 
dense environments is seen, at z = 1 the density-morphology relation becomes fiatter 
and late type galaxies are even seen to have a higher clustering amplitude than early 
types. 
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1 INTRODUCTION 

Forming realistic galaxies is one of the grand challenges for 
numerical cosmology. A promising approach assumes that 
this can be done, given the correct cosmological model and 
enough spatial and mass resolution in a simulation of the 
physics of gravitational collapse of gas and dark matter. Pi- 
oneering work in this vein (e.g., Evrard 1988, Hernquist & 
Katz 1989) included some of the physical implementations 
(Smoothed Particle Hydrodynamics, Gingold & Monaghan, 
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1977, Lucy 1977) still most used today. Gradual improve- 
ments and additions, including star formation (e.g. Katz 
1992), feedback (e.g., Springel and Hernquist 2003a), black 
hole accretion (e.g., Di Matteo et al. 2005, Springel et al. 
2005, Okamoto et al. 2007) have led to better modelling of 
galaxy properties (e.g., Abadi et al. 2003a,b, Governato et al. 
2007). The most striking result to come from this research 
effort is the generic problem of too little angular momen- 
tum and too centrally concentrated galaxies (see e.g., the 
recent review by Mayer et al. 2008). Many aspects of our 
understanding of galaxy formation have been progressing, 
however, spurred on by high resolution "zoomed" resimula- 
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tions of individual galaxies (e.g. Abadi et al. 2003, Sommor- 
Larsen ct al. 2003, Governato et al. 2004, Robertson et al. . 
2004, Zavala et al. 2007) In the meantime, larger-scale sim- 
ulations (e.g., Springel et al. 2005, Weinberg et al. 2004) 
have led to better predictions of the large scale clustering 
properties of galaxies. In the present paper, we aim to ad- 
dress questions which link both approaches, by looking at 
internal properties of galaxies which form in a cosmological 
volume with uniformly high resolution. In this way we can 
look at not only the properties of a wide range of galaxies, 
but the relation of these properties with their large scale 
environments. 

Using the same simulation to study individual galax- 
ies and the evolution of structure in a cosmological volume 
is necessary to ensure that we deal with a fair sample of 
galaxies. In zoomed resimulations (e.g., Katz et al. 1994, 
Robertson et al. 2005) galaxies that are chosen to be followed 
at high resolution (usually systems relatively free of major 
mergers at late types) may be special by virtue of their se- 
lection, and it is not easy to quantify biases that this may in- 
duce in statistical studies of their properties. Running a cos- 
mological volume with enough resolution to carry this out to 
rcdshift z — is extremely expensive computationally. Mo- 
tivated both by this and the rapid growth of observational 
data available for redshift « = 1, (e.g., the DEEP2 study of 
galaxy properties and environment at « ~ 1, Cooper et al. 
2007, the VVDS survey, Meneux et al. 2008, and COMBO- 
17 survey. Bell et al. 2004) we focus most of our attention 
on a large, high resolution simulation run to this redshift. 
Results from this computation, known as the BHCosmo run 
were first presented by Di Mattco ct al. (2007). It features 
the first direct inclusion of supermassivc black hole accretion 
and feedback in a cosmological hydrodynamic model. 

The questions which can be asked include whether we 
are able to make disk galaxies which look like observed spi- 
rals, and whether they lie on the TuUy-Fisher relation. Hav- 
ing a large statistical sample of simulated galaxies will en- 
able studies of how their properties are correlated and in- 
terrelated, for example how disk scales lengths are related 
to circular velocities and so on. In the present study we do 
not carry out population synthesis modelling, instead con- 
centrating on morphology of stellar mass density, stellar and 
gas kinematics and other related measures of galactic struc- 
ture. 

Obscrvationally, large datasets of galaxies at redshifts 
z ^ 1 arc becoming available. As mentioned above, with the 
DEEP2 redshift survey (see e.g.. Coil et al. 2008 and ref- 
erences therein) data it has become possible to look at the 
dependence of galaxy clustering on color, luminosity and 
other properties. HST images of galaxies with enough reso- 
lution for morphological classification have extended studies 
of the morphology-density relation to z ^ 0.5 — 1.0 (see 
e.g., the early work of the MORPHS collaboration. Dressier 
et al. 1997 through Smith ct al. 2005). At higher redshifts, 
Peter et al. (2007) have also measured galaxy morphologies 
in the environment of a rich cluster at redshift 2: = 2.3 us- 
ing HST/ACS images. With the increase in observational 
data available at these redshifts, it is important to inves- 
tigate the robustness of theoretical modelling with a vari- 
ety of techniques. For example, Zheng et al. (2007) have fit 
Halo Occupation Distribution models to the luminosity de- 
pendence of clustering in DEEP2 and Croton et al. (2006) 



have made semi-analytic model predictions for the cluster- 
ing strengths of red and blue galaxies at redshifts relevant 
for DEEP2. In the present work, we explore the predictions 
of hydrodynamic simulations. 

Our plan for this paper is as follows. In Section 2 we 
describe the cosmological simulations, how we select galax- 
ies from the simulation outputs and show some images of 
examples of late and early type galaxies. In Section 3 we 
investigate galaxy properties in detail, including kinemati- 
cally defined bulge to disk ratios, and projected stellar den- 
sity profiles. We show how measurements of these quantities 
are affected by simulation resolution and then examine the 
relationships between all the properties for individual galax- 
ies. We also plot the angular momentum content of late type 
galaxies and compare it to observations, as well as compar- 
ing to the observed stellar mass TuUy-Fisher relation. In Sec- 
tion 4 we focus on the relationship between galaxy properties 
and environment, including comparison to the observed den- 
sity morphology relation. In Section 5 we look at large scale 
structure in the simulated galaxy distribution and compute 
the autocorrelation function. We summarize and discuss our 
findings in Section 6. 



2 SIMULATION 

We have carried out a large cosmological simulation, using 
the SPH code Gadget-2 (Springel 2005), with the inclusion 
of modelling of black hole accretion and feedback, previously 
used in simulations of isolated galaxies and galaxy mergers 
(Di Mattco ct al. 2005, Springel et al. 2005). This direct cos- 
mological simulation of black hole evolution in a cosmologi- 
cal context, known as the BHCosmo simulation is presented 
and described in more detail in Di Matteo et al. (2007). 
The simulation algorithms and description of the black hole 
population that forms and its relationship to galaxies in the 
simulation are also set out in that paper. A slightly different 
approach to modelling black holes in cosmological Gadget-2 
simulations, including feedback from a "radio mode" is also 
presented by Sijacki et al. (2007). Okamoto et al. (2007) have 
also incorporated their own modelling of black hole physics 
in zoomed simulations (also run with Gadget-2) of individ- 
ual galaxies. We discuss some similarities with our results in 
Section 6.2 

Besides black hole feedback, the simulation includes the 
multiphase star formation model of Springel & Hernquist 
(2003a), in which stellar wind particles can propagate from 
the sites of star formation and carry energy over large dis- 
tances. This model was previously used in the zoomed simu- 
lation of an individual disk galaxy in a cosmological context 
by Robertson et al. (2005). 

The cosmological parameters used as inputs to the BH- 
Cosmo simulation are consistent with the WMAP first year 
results (Spergel et al. 2003), i.e. h = 0.7 =H()/100 km s"^ 
Mpc~^, matter content Qrn = 0.3, cosmological constant 
flA = 0.7, baryon content fit = 0.04 and amplitude of mass 
fluctuations, as = 0.9. The simulation volume is a periodic 
cube of side length 33.75 ft~^Mpc, and the BHCosmo con- 
tains 486'' dark matter and initially 486^ gas particles, giving 
a mass resolution of 2.75 x lO'^h~^M0 per dark matter par- 
ticle and 4.24 x lO®h~^M0 per gas particle. We evolve the 
simulation from redshift z = 99 down to « = 1, using con- 
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stant gravitational softening in comoving units. At z = 1, 
the force resolution is 1.4 proper /i~^kpc. 

The BHCosmo simulation has a mass and force resolu- 
tion that fits into the scheme devised by Springel & Hern- 
quist (2003b), who ran a set of simulations with box size 
varying by factors of ~ 3 and particle numbers by factors of 
1.5^. In their naming scheme, the BHCosmo would be equiv- 
alent to a D6 simulation. We have also run another simula- 
tion with the same initial conditions, box size and physical 
modelling but coarser mass and force resolution (equivalent 
to D4, i.e. 2 x 216^ particles). We use this D4 run, which has 
a mass per particle 11.4 times greater than the BHcosmo, 
and force resolution 2.3 times larger, to gauge the effects of 
resolution on our results, as was done by Di Matteo et al. 
(2007). Because the D4 was run with the same random initial 
Fourier phases as the BHcosmo it is possible to identify the 
same objects which form in both runs and therefore make 
direct comparisons (we shall do this in §3.3). 

2.1 Galaxy selection 

We select galaxies from the particle distributions using a 
variant of the subgroup finder SUBFIND (see Springel et al. 
2001). This selects gravitationally bound clumps of particles, 
using the gas, dark matter and stellar content of groups to 
do this. Each subgroup is therefore deemed to be an individ- 
ual galaxy. Each galaxy is also associated with a friends-of- 
friends group, which is found using the usual method, link- 
ing together particles connected by less than 0.2 times the 
mean interparticle separation. Particles within the friends- 
of-friends group are approximately those that lie within the 
virial overdensity threshold of a dark matter halo. 

From the list of galaxies, we pick those which have more 
than 5000 particles in the BHCosmo run, in order to limit 
the noise on morphological and kinematic properties mea- 
sured from the galaxy particle distributions. At z = 1, we 
have 1180 galaxies in the BHCosmo simulation. The parti- 
cle number cut is essentially equivalent to a cut in circular 
velocity of 100 km s~^. We find that aX z = 1, galaxies of 
similar mass to the Milky Way (those with circular velocities 
within 10% of 220 km s"^) contain ~ 60000 particles on av- 
erage. For the D4 low resolution run we use an appropriately 
lower threshold on the particle number (439 particles). 

For each galaxy (i.e. subgroup) we define a central 
point, which is the particle (star, gas or dark matter) at 
the minimum of the gravitational potential. This particle 
serves as a center when computing the stellar density and 
angular momentum profiles and for all other purposes in the 
paper when the center of the galaxy is required. 

2.2 Example galaxies 

Before characterizing the statistical properties of our sample 
of galaxies, we examine their morphology visually. To do 
this, we take galaxies from the redshift z = 1 output of the 
BHCosmo run. In Figures [T] and [2] we show 10 representative 
galaxies taken from two classes, those with high D/T (Disk 
to Total) ratios and low D/T ratios respectively. The D/T 
ratio is derived from a kinematic decomposition of the stellar 
mass content of each galaxy (following e.g., Abadi et al. 
2003a), and will be more fully explored in §3.1 For now it is 



sufficient to state that these decompositions form the basis 
for one possible separation into "late" type disky galaxies 
and "early" type ellipticals. 

Figure [T] shows galaxies that lie principally in the top 
10% or so of D/T values, although because most galaxies are 
bulge dominated, even these include some with D/T values 
less than 0.5. We note that zoomed resimulations of disk 
galaxies in other work typically have similar D/T ratios (see 
e.g., Abadi et al. 2003b in which the galaxy had D/B of 
0.83, or D/T of 0.45). The plot shows two views each of the 
stellar density, gas density and dark matter density, on a log- 
arithmic scale. If the angular momentum vector of the stellar 
component defines the z-axis, the two panels for each compo- 
nent show the face on [x — y) view and a side on {x — z) view. 
We can see that the stellar components are fiattened in the 
direction expected if they are rotationally supported. They 
are not particularly thin, and an obvious bulge is not easy 
to spot in the surface brightness plots. Comparing with the 
zoomed simulations of disk galaxies performed by Robertson 
et al. 2004 (their figure 3, top left) or Governato et al. 2007 
(their figure 6, left panels), we can see similar morphologies. 

The gas components in Figure [T] are more varied in ap- 
pearance than the stars, with most galaxies having relatively 
thin disks. Some however (such as the second and third from 
bottom) have rather clumpy distributions and are less ob- 
viously flat. Approximately half the dark matter halos can 
be seen to be obviously slightly flattened in the same direc- 
tion as the stellar disks. Alongside the plots, we have listed 
the D/T ratios and the total mass of the subhalo. Most of 
the disk galaxies we plot are fairly large, with total mass 
around 10^^h~^MQ. The fifth from the top, the largest and 
most massive has the lowest D/T ratio. We also show next 
to each panel an indicator of the physical size of the galaxy, 
the radius of a sphere which encloses half the stellar mass 
(riiaif), in proper ft^^kpc. 

Moving to the early types in Figure [2] one immediate 
difference that can be seen from the stellar density plots is 
that they do appear to be more centrally concentrated than 
the disks. On average, the ratio of major to minor axes ap- 
pears to be lower, and the orientation of the galaxy also 
appears to have little relation to the angular momentum di- 
rection (as expected, give that these galaxies were selected 
on the basis of their stellar kinematics). The gas compo- 
nents are almost all rather clumpy and roughly spherically 
distributed, although the first galaxy plotted has a very ob- 
vious thin gas disk. This is rather unusual, given the almost 
spherical nature of its stellar distribution, although we note 
that it does have the highest D/T ratio of all the galaxies 
in this plot. Interestingly, the dark matter contours do not 
appear to be oriented in the same way as the stellar dis- 
tribution, as any visible degree of correlation between the 
major axes of the dark matter and stars is hard to pick 
up on. The gas distributions do appear perhaps to line up 
slightly better with the dark matter. We have however com- 
puted the median angular difference between the directions 
of the principal axes of the stars and gas (we return to this 
in §3.7 below) and find a value of 24 °, a fairly strong align- 
ment. This value is for all galaxies though, and presumably 
strongly influenced by the good alignments of gas and stellar 
disks. 
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Figure 1. Images of 2 = 1 galaxies from the BHCosmo simulation selected kinematically to have a high D/T ratio (i.e., "late types", 
see text §3.1). We show 10 galaxies, in no particular order, from the top to the bottom of the page. The three columns with different 
appearances are from left to right the stellar surface density, the gas surface density and the dark matter surface density in a cube 
centered on the minimum of the gravitational potential for the galaxy (see §2.1). In each column, there are two panels corresponding to 
two projections, the leftmost is such that the galaxy is face on, using the angular momentum vector of the star particles in the galaxy 
to define the direction. The rightmost panel is an edge on projection. All surface densities are plotted with a logarithmic scale, and the 
contours in the dark matter plots are shown at intervals of 0.2 dex The length scales are such that each panel shows a region of space 
from — 2riiaif to +2riiaif, where rjj^if is the radius which encloses half the stellar mass. The values of rj^^if for each galaxy (in proper 
^"'^kpc) are shown to the right, along with the D/T ratio and total mass of the galaxy /subgroup (including stars, gas and dark matter.) 



3 GALAXY PROPERTIES 

We will focus on three main properties of our simulated 
galaxy sample, their kinematics (angular momentum and 
kinematic decomposition into disk and bulge), their profiles 
(quantified by a Sersic index) and their mean stellar ages and 
star formation rates. We do not classify their morphology by 
eye as has been done observationally (e.g., by the MORPHS 



group. Dressier et al. 1997 ), but we will use the kinematic 
decomposition and Sersic indices as indicators of morphol- 
ogy, for example when investigating the morphology-density 
relation. We make use of our two simulations, the BHCosmo 
run and the D4 run which have identical initial conditions 
but different mass and spatial resolution to investigate the 
effects of biases in galaxy properties due to numerical effects. 
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Figure 2. Images of ten 2 = 1 galaxies from the BHCosmo simulation selected kinematically to have a low D/T ratio (i.e., "early 
types", see text §3.1). We show the stellar, gas and dark matter surface density from left to right, in the same fashion as the late type 
galaxies in Figure [T] 



Comparisons of both simulations are presented in detail in 
§3.3 below, including some quantification of the effects of 
resolution. We first explain the measures of kinematics and 
stellar profiles, also including examples drawn from both 
simulations. 



3.1 Kinematics and D/T ratios 

We examine the stellar orbits in the galaxies, following a 
technique often used in the literature (e.g., Abadi et al. 
2003b), we use the angular momenta of individual stars 
to dynamically decompose the galaxies into disks and 
spheroids. 

In the disk galaxies (Figure [1]), the total angular mo- 



mentum vector of the stars defines a direction which is vi- 
sually perpendicular to the stellar disk. We will investigate 
how well this holds quantitatively later (§3.7), but for now 
we use this direction to define the z— axis. For a star par- 
ticle, we compute the specific angular momentum in the z 
direction, J and also its specific binding energy E. Following 
exactly the procedure of Abadi et al. (2003b), we compute 
E relative to the total mass of the galaxy within the viral 
radius, r2oo . For each star particle we then compute the spe- 
cific angular momentum of the corotating circular orbit with 
the same binding energy, Jc- The ratio J/Jc gives for each 
particle a measure of the orbital circularity. The maximum 
value, J/Jc = 1 occurs for particles which are following in 
circular orbits in the same direction as the majority of the 
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Figure 3. Histograms of the orbital circularity of star particles, 
J/ Jq for the galaxies (at 2 = 1) shown in Figures [T] (late types, 
left column) and|2](early types, right column). The quantity J/Jc 
is the ratio of the orbital angular momentum of an individual star 
particle to one with the same energy but on a circular orbit (see 
text). The shaded histograms are for galaxies selected from the 
BHCosmo simulation and the superimposed lines are the same 
galaxies taken from the lower resolution simulation, D4. How the 
matching galaxies were found in the D4 run is described in §3.2. 
The y-axis in each case is the relative probability density, nor- 
malized to 1 at the maximum value. 



stars. Particles with close to this value are therefore likely 
to be in a disk component. 

In Figure O we show histograms of the J/ Jc values for 
all stars for our 20 example galaxies, with the late types on 
the left, and early types on the right. In order to avoid going 
out very far into the outskirts of the galaxy where particle 
numbers are sparse, we only plot particles which lie within 4 
half stellar mass radii (rhaif) of the center of the galaxy (the 
results are insensitive to this choice). To compute the D/T 
ratio, we make the (somewhat arbitrary) assumption that 
the bulge has zero total angular momentum (see e.g., van 
den Bosch et al. 2002). Starting from the negative end of 



the J/ Jc distribution we associate all the star particles with 
negative J/Jc with the bulge as well as an equal mass of 
star particles with the equivalent positive J/Jc values. The 
star particles left over then make up the disk. In this way, 
we compute the disk to total ratio, D/T already enumerated 
next to the galaxy plots in Figures [T] and [5] 

Because the bulge/disk decomposition relies on this 
J/Jc distribution it is not surprising that the galaxies in 
the left panel, chosen to be more disky have asymmetric his- 
tograms. Eight of the ten disk galaxies have the highest point 
in the J/Jc distribution associated with the region around 
J/Jc ~ 1, showing that the galaxies indeed have a strong 
stellar rotationally supported component, as we would ex- 
pect from their morphology. Some of the distributions are 
markedly more asymetric than others, with 2 not showing 
any real evidence for a spheroid component (a bump cen- 
tered on J/Jc ~ 0). The shaded histogram is the result for 
our high resolution BHCosmo simulation, and the (some- 
what noisier) line is the result for the same galaxies taken 
from the low resolution D4 model (see Section 3.3 for details 
on the matching of galaxies in the simulations). We can see 
that the two simulations can give rather different results for 
the J/Jc distribution for the same galaxy. This is likely to 
be due among other things to stochasticity coming from ini- 
tial fluctuation power below the Nyquist limit of the lower 
resolution model, and differences in the numerical integra- 
tion of particle orbits coming from resolution effects. This 
will lead to a scatter in D/T values between the two simula- 
tions. We investigate quantitatively the effect of resolution 
on this parameter below. 

The low D/T galaxies in the right panel of Figure |3] un- 
surprisingly have central peaks in the J/ Jc distribution. The 
effects of resolution on the histograms themselves appears to 
be somewhat lower, and in fact we shall see that in fact the 
scatter in D/T values from the two simulations is closer to 
being constant in log D/T. Interestingly, if we search for ev- 
idence of a small stellar disk to go along with the gas disk 
seen in the visual inspection of the first galaxy, we see that 
there is no evidence for one in the J/Jc distribution. 

Although we do not plot J/Jc for the gas particles, we 
have examined their distribution and find it to be more ro- 
tationally supported than the stars, in agreement with pre- 
vious work (see e.g., Mayer et al. 2008). The previous work 
on resimulated spiral galaxies has yielded D/T values from 
kinematic decomposition of e.g., 0.41 (Abadi et al. 2003b), 
0.67-0.76 (Okamoto et al. 2007) 0.74 (Governato et al. 2004). 
Our example disk galaxies shown here are generally consis- 
tent with this range. From our large sample, we can say 
something about how rare truly disk dominated galaxies are 
in the simulation. At z = 1 only 3.5% of the galaxies in the 
BHCosmo run have D/T greater than 0.5. The maximum 
D/T value we find is 0.76 (a D/B ratio of 3.22). This is the 
same as the (mass weighted) D/T ratio of the disk galaxy 
simulated by Okamoto et al. (2007), which included AGN 
feedback. 

Observationally, spiral galaxies in the local Universe are 
seen to have D/T ratios from unity (bulgeless disks) to ~ 0.1 
(see for example Graham et al. 2001, who find a median 
D/T value of 0.75 for a sample of 86 face-on spirals). The 
D/T ratios are usually (and in the case of Graham et al. ) 
calculated in terms of luminosity of each component, using 
fits to the photometry, rather than in stellar mass from the 
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kinematics, as we have done. This will tend to make our D/T 
rather lower because of the relative mass/light ratio of disks. 
Nevertheless, it is clear that the disk components of our 
galaxies are small, a fact related to the angular momentum 
problem in simulations. In separating galaxies into early and 
late types, we will therefore use a rather low value of D/T 
to divide the two classes (we use D/T=0.2.) 

3.2 Projected stellar density profiles 

There are a number of ways to decide on galaxy type 
using fits to the radially averaged luminosity profiles of 
galaxies. Galaxy classification based on photometry often 
makes use of the fact that early type galaxies are more cen- 
trally concentrated. For example, elliptical galaxies and the 
bulges of spirals tend to obey the deVaucouleurs profile well: 
E oc exp (r/ro)"^'^^, and disk components have pure expo- 
nential profiles, E oc exp— r/rcff (where E is the projected 
stellar density). A common method for using this informa- 
tion is to fit a Sersic profile which can account for both, 

E = exp(r/ro)-'''" (1) 

where n is the Sersic index (Sersic, 1968). At redshift z = 
the largest study of Sersic index fitting was carried out by 
Blanton et al. (2005). At z ~ 1, there are quantitative mea- 
surements of the morphology from a variety of surveys in- 
cluding the CFHT legacy survey (Nuijten et al. 2005,, using 
the Sersic index) COSMOS (e.g., Capak et al. 2007, using 
the Gini parameter). 

To examine the stellar profiles of our simulated galax- 
ies, we rotate them so that they are face on with respect to 
the stellar angular momentum vector. We then use this pro- 
jected distribution to bin the stellar particles, assigning their 
masses to radial bins. Just as with the kinematics above, we 
are using stellar density rather than luminosity in a partic- 
ular band. 

In Figure|4]we show the projected stellar density profiles 
for our 20 example galaxies, again in the same order from 
top to bottom and with the late types on the left. We show 
results for the D4 run as a dashed line as well as the BH- 
Cosmo run in solid. The differences between the profiles for 
the two simulations appear to be much smaller than those 
in the kinematic decompositions in Figure |3l However, we 
shall see later that these small differences actually lead to 
large scatter in the best fit Sersic index for each galaxy in 
the two simulations. 

In the left panels of Figure |4] we have also plotted an 
exponential disk profile, fitting it to the radial bins at greater 
distances than rhaif from the center of the galaxy. We ignore 
the central bulge region within r^aif in order to void biasing 
the fit, and also to avoid fitting the central profile inside the 
gravitational softening length of the simulation (1.4 proper 
h~^kpc for the BHCosmo). Looking at how well the late 
type galaxies obey a pure exponential profile, we can see 
that 4 out of the 10 have no significant bulge component, 
with the 6th from top being the best. A galaxy which is 
pure exponential all the way to the center was seen in the 
resimulation of Robertson et al. (2004) . Here we see that this 
can happen relatively often, even with the comparatively low 
resolution of the D4 simulation. 

In most of the panels in Figured we can see that in the 
inner regions, the D4 simulation falls somewhat below the 
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Figure 4. Projected stellar density profiles of the galaxies shown 
in Figure [T] (left column, late types) and Figure [2l (right column, 
early types) . The thin solid lines are the profiles from the galaxies 
in the BHCosmo simulation and the dashed lines are the profiles 
for the same galaxies taken from the lower resolution D4 simu- 
lation. The units for E, the projected stellar density are h~^MQ 
(proper kpc)~^. The x-axis is radial distance from the center of 
the galaxy, as a fraction of the half stellar mass radius (given for 
each galaxy in Figures [T] and [2] In the left panels we also show as 
thick grey lines exponential profile fits to the profiles (beyond the 
^'half radius) of the late type galaxies. In the right hand panels 
we show, also as thick grey lines de Vaucouleurs profile fits to the 
early-type galaxies. 

higher resolution run. As found for example in the conver- 
gence tests of Governato et al. (2007), the inner parts of the 
galaxies are not resolved, so that we cannot say much about 
the structure of the bulges. Governato et al. (2007) esti- 
mate that several million resolution elements will be needed 
within the virial radius to do this. 

Looking at the early types in the right panel, we can 
see that they are on average more centrally concentrated, as 
we noticed from looking at the images (Figure [2]) We have 
plotted a de Vaucouleurs profile fit to the profiles in these 
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images, and we can see that for 9 out of the 10 early types 
it gives reasonable results. We also see however that one of 
them (the 3rd from bottom) is also pretty close to exponen- 
tial in shape, and the de Vaucouleurs profile is not a good 
fit at all. Returning to the image of this galaxy in Figure [2j 
it can perhaps be seen to be fairly smooth in the center, but 
there is no real hint that its profile is peculiar, nor is there in 
the kinematics Figure |3}- This is quite suprising, and shows 
that even a good exponential profile is not necessarily a sign 
of a disk dominated galaxy. 

In section 3.4, we shall instead fit a Sersic profile to all 
the galaxies, and use that to make a separation of galaxies 
into two classes. For now we have seen that for most of our 20 
example galaxies there is reasonable agreement between the 
type suggested by their kinematics and their stellar profiles. 



3.3 Resolution tests 

We have mentioned two of the galaxy properties we are most 
interested in above, the D/T ratio and the Sersic index, n. 
We will also study the circular velocity, Kirc, the total an- 
gular momentum of different galaxy components (e.g., jgas), 
the star formation rate, mean star formation redshift (zform) 
and the exponential disk scale length, rdiak- For all these 
quantities, the numerical resolution of the simulation could 
infiuence the results. In this section we investigate whether 
changing the resolution causes a systematic bias by compar- 
ing the same galaxies in the BHcosmo and D4 runs. Because 
these simulations were run with the same intitial random 
seeds to set up their initial conditions and have the same 
box size, in principle we should be able to find high and 
low resolution versions of the same galaxies. With such a 
resolution test, we will be able to get information on which 
of the parameters is most reliably computed by the simula- 
tion, in terms of both systematic bias with resolution and 
amount of scatter in results between the two resolutions. In 
order to emphasise the resolution difference between the two 
simulations, we refer to the BHcosmo run as the D6 run in 
this section (using the nomenclature of Springel & Hernquist 
2003). 

In order to carry out this test we first need to match 
up the galaxies in the D4 and D6 runs. We first apply a 
particle number cut which is equivalent to the 5000 particle 
lower limit for the D6 run. After applying this cut (439) 
particles we are left with 1297 D4 galaxies. For each D6 
galaxy, we look for the closest (the smallest 3 dimensional 
separation between the centers, defined by the particle at 
the minimum potential) one in the D4 output. We also apply 
another criterion, so that the galaxies must be within 25 % 
of each other in total mass to be considered a match. We do 
not exclude galaxies which have already been matched. The 
procedure is insensitive to this, as is it is to changing the 
25% matching criterion between 10% and 30%. In order to 
gauge roughly how many galaxies have been mismatched, we 
try applying the same mass similarity criterion to the dark 
matter masses of the galaxies. We find that none of the top 
20 galaxies by mass change, but below that approximately 
10% of the D6 galaxies are matched to another in the D4. 
We find similar results with the stellar mass, and so our 
conclusion is that approximately 10% of the galaxies are 
likely to be mismatched. This is a small number, but should 



be borne in mind when considering the comparisons between 
simulations. 

In Figure O we show the results for 9 different quanti- 
ties, plotting the value for the D4 simulation against that 
in the D6 simulation, for each of the 1180 galaxies. In the 
plots we show two categories of galaxy, those which have 
> 439 particles in the D4 (equivalent to 5000 in the D6) as 
open symbols and those which have > 5000 particles in the 
D4 (equivalent to 57000 in the D6) as filled points. This is 
so that we can see whether any scatter or bias is less for 
larger particle numbers per galaxy. As the D6 galaxies all 
have > 5000 particles, in this sense the filled points will give 
a sense of the uncertainties induced by resolution on the 
results in the rest of the paper, which use the D6 outputs. 

Looking at the different panels of Figure [S] we can see 
that all the galaxy parameters show correlations between 
the two simulations, although the scatter between them is 
in some cases (e.g., Sersic index, n) much larger than oth- 
ers, (e.g., V^irc). An X = y line has been drawn on all panels 
and the first impression is that while the scatter may be 
large, there is not much evidence of strong systematic bias 
caused by decreased resolution. The most obviously differ- 
ent slope from 1 is that in the Zform panel. Here we can see 
that galaxies in the D4 simulation tend to form their stars 
on average significantly later than in the D6. For example, 
stars which formed at average z = 3.5 in the higher res- 
olution model form at average z = 2.5 in the D4 run (for 
galaxies with between 439 and 5000 particles in the D4) . The 
effect of the bias is seen to be significantly reduced when we 
move to galaxies with higher numbers of particles (the filled 
points). This bias with resolution is expected, given that 
at higher resolution smaller scale density fiuctuations are 
present which can collapse earlier. A detailed study of star 
formation in the same multiphase model used in the present 
simulation was carried out by Springel & Hernquist (2003a) 
and to which we refer the reader for further resolution tests 
and interpretation. 

Related to the star formation history is the current SFR 
at 2; = 1, which is also shown in Figure [S] in this case we 
can see that the SFR in the D4 run is systematically higher 
(the points tend to lie above the x = y line) than for the 
D6 galaxies. This is particularly true for the low particle 
number galaxies. A possible interpretation is that perturba- 
tions which were not present in the D4 compared to D6 due 
to insufficient resolution are overtaken by collapse on larger 
scales and forming stars. 

The two parameters which have the most scatter are the 
D/T ratio and the Sersic index n. The strength of the corre- 
lation between the two simulations can be judged using the 
Pearson correlation coefficient, which for D/T is 0.22 for all 
galaxies and 0.35 for D4 galaxies with > 5000 particles. The 
significance levels for the correlations are 10"^* and 2x10"^ 
respectively (there are ~ 10 times fewer galaxies for the high 
particle number subset, so the significance is lower). If we 
fit a straight line y = a + bx to the D/T results we find an 
offset of a = 0.21,0.11 for the small and large galaxies, re- 
spectively. This means that there is a small systematic offset, 
but it is in the direction of the low resolution galaxies having 
higher D/T ratios, not what one would expect if simple lack 
of resolution causes the angular momentum problem. This 
fit assumes that all the errors lie in the coordinate (the 
D4 D/T ratio) and gives slopes of 6 = 0.54 and b = 0.60. 
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Figure 5. Comparison of galaxy properties in the BHcosmo (referred to here as the D6) simulation and the same galaxies selected from 
the lower resolution D4 run (see §3.3 for details of how the matching galaxies were found). In each panel, the open circles show galaxies 
which contain > 5000 particles in the D6 run, which corresponds to an equivalent of > 439 particles in the D4 run. The filled points 
are for galaxies which contain > 5000 particles in the D4 run (or the equivalent of 57000 for the D6). In the top row, we show from left 
to right the circular velocities of galaxies (measured at the 2.2 times exponential scale radius rjjgj- of the galaxies), the kinematically 
determined D/T ratio, the mean formation redshift of star particles in the galaxy. In the middle row we show the star formation rate, 
the Sersic index, and the exponential scale radius. In the bottom row are plotted the mean specific angular momenta of the stars, the 
gas and the dark matter. 



The fact that the slope is not equal to one may be partly 
due to a significant portion of tlie scatter not coming just 
from the y— coordinate (i.e., the D6 value is inaccurate) and 
partly due to resolution lowering the angular momentum. 
Because the relevant "errors" in the x and y-coordinates are 
not known, it is not possible to say more than this. 

For the Sersic index, we have a significantly worse cor- 
relation for all galaxies, r = 0.11, with a significance of 
8 X 10~^. For the > 5000 particle galaxies it is surpris- 
ingly good, r = 0.84 (significance < 10"^""). Visually, the 
scatter does appear to be very large, something that is in- 
teresting given the good correspondence between the lines 
showing the stellar profiles in the D4 and D6 cases in our 
plots of example galaxies (Figure |4]). This shows that dif- 
ferences between the slopes of the profiles are quite subtle 
and that n may be a problematic indicator of galaxy type. 



Given the evidence of the profiles seen in Figure |4j the phys- 
ical gravitional softening length difference between the two 
simulations, (1.4 ft~^kpc vs 3.1 h~^kpc) does not seem to 
make a difference to the Sersic indices. 

The bottom three panels of Figure [S] show the reso- 
lution comparison of the specific angular momentum j of 
each of the three components, gas, dark matter and stars, 
of each galaxy. The scatter is relatively large, worse visually 
than the impression produced by K:irc and by rdisk- It does 
decrease with higher particle number and there does not ap- 
pear to be a significant bias with resolution. This may mean 
that increasing resolution is not the whole story for dealing 
with the angular momentum (see e.g., Agertz et al. 2007) 
or else that with much higher resolution we might suddenly 
see a dramatic change (see e.g., Ceverino & Klypin 2007 for 
exploration along these lines) . One fact which is obvious is 
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Figure 6. Galaxy properties in the BHCosmo simulation at redshift 2 = 1. In the panels which contain scatter plots, we show 1 point 
for each galaxy. We also show the median value in bins as a thick solid line, with Poisson error bars. The 90th and 10th percentile 
distibutions in the same bins as the median are also shown as thin solid lines. The thick dashed line in each case shows results for the 
galaxies taken from the lower resolution D4 simulation. The diagonal panels show the probability distribution of the parameter values, 
on the same x-axis scale as the scatter plots. Again the dashed lines show the pdf from the D4 run and the solid line enclosing a shaded 
region the results from the BHCosmo simulation. The y-axis of the diagonal panels represents a logarithmic scale of probability. The 
overall normalization is arbitary, but the height of each panel is 4 dex. 



that the gas has higher j values than either the stars or the 
dark matter. We shall return to this in Section 3.5 below. 

The conclusion we can draw from our study of identical 
galaxies at lower resolution is that most of the properties are 
robust, although the scatter between two resultions can be 
large. It is largest for the Sersic index n and D/T ratios, two 
of the parameters we will use to separate our galaxies into 
different types. A relatively large systematic bias between 
parameters due to resolution effects is only present for mean 
stellar age, which is affected in the way one expects due to 
missing small scale power. 



3.4 Correlations between galaxy properties 



One of the aims of this paper is to see whether galaxy mor- 
phology in simulations is affected by environment. The other 
main aim is to see how galaxy properties are interrelated. 
For example, we look at whether disk dominated galaxies 
are smaller than bulge dominated ones, how Sersic index is 
related to D/T ratio and so on. Our main tool in this regard 
is a scatter plot of 6 important galaxy properties where each 
is plotted against the 5 others. This is shown, for 2 = 1, the 
redshift we are focussing on (for the BHCosmo run) in Fig- 
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ure[6] For each of the quantities we show a histogram of the 
values on the diagonal (where the y— axis is in log units). 

We have seen in Figure |4] that it is possible to roughly 
fit an exponential profile to the region outside the inner half 
stellar mass radius, even when the galaxy is an early type. 
In this way, we have computed an effective radius r^s for all 
galaxies, so that we can compare early and late types on the 
same footing. We use this value in our measurement of the 
galaxy circular velocity, which we compute using 

Vcirc = A/GM<2.2r„ff/roff. (2) 

Here M<2.2r„ff is the total mass within 2.2 effective radii 
(e.g., Courteau 1997). Our Kirc values for all galaxies are 
therefore computed in this way. We use reff in order to be 
consistent with observational definitions of circular velocity, 
but find no significant difference to our results if we instead 
use rhaif. 

In the panels of Figure [S] we show lines which represent 
the medians of the points in bins, as well as the 10th and 
90th percentile lines. For the median, as well as the pdfs, we 
show results for the D4 run as a dashed line. We note that 
half of the panels are repeated, in a different orientation. 
We include them because seeing the median of the other 
quantity is still useful. 

A number of interesting facts are apparent from a look 
at Figure [6] The narrowest relation is between reff and V^irc. 
Galaxies with Vcirc = 200 km s~^ have a median r^s of 4 
proper /i~^kpc, or 5.7 kpc. This can be compared to the 
Milky way, which has reff = 3.5 kpc, for example. It seems 
that the size of galaxies in the simulation is therefore roughly 
correct, at least within a factor of 2. The main problem, as 
we shall see from the D/T ratios is that they have bulges 
which are too large. The values of reff and Vdic can be seen 
from Figure [6] to be unaffected by resolution. 

Moving on to the D/T ratios, we can see that there 
is a trend for disk dominated galaxies to be smaller (lower 
Kirc and reff values) than systems with greater bulge compo- 
nents. The median D/T value for the smallest galaxies plot- 
ted (Vcirc ~ 100 km s~^) is approximately twice that of the 
largest {Vdrc ~ 500 km s^^). This trend can also be seen in 
the 90% envelope, and the 10% envelope, which both trend 
downwards with galaxy size. There is a particularly obvi- 
ous lack of disk galaxies above Virc ~ 300 km s~^. Even at 
z — 1, therefore the simulations exhibit the well known ob- 
servational trend for the largest galaxies to be early types. 
While the trend of D/T with size does seem to be reasonable 
compared to observations, the absolute values are low, as we 
have seen before. For example, there is a lack of disk galax- 
ies which have no detectable bulge, but which are known 
to exist in the local Universe (e.g., Dalcanton & Bernstein 
2000). 

If we look at the age of galaxies next, it becomes ap- 
parent that the behaviour of galaxy property trends in the 
presence of large scatter is quite complex. For example in 
the panel of Virc vs age, the higher Vcirc galaxies are slightly 
older in the median, although the errors are large. On the 
other hand, in r^s vs age, the smaller galaxies (these with 
Teff < 5/i~^kpc) have higher ages. Conversely, the lower D/T 
galaxies also have higher ages, as we would expect of late 
types. We find therefore that some of the usually assumed 
relationships between galaxy properties and type (see e.g.. 
Mailer et al. 2006, Koda et al. 2007) hold, but in the grid 



there are several counterexamples. The effect of resolution 
on galaxy ages is for them to be lower in the D4 run, as we 
have seen in Section 3.3. This effect is shown by the median 
line in the panels of Figure |6l We do see however that al- 
though the ages are lowered by a constant offset, the trend 
of the ages with other properties is broadly similar to the 
higher resoution run in all the panels. 

The Sersic index, n, of the galaxies can be seen to in- 
crease with reff (although the last r^g bins drops down, the 
errors are large). Diskier galaxies also have slightly lower n 
values, as we would expect, with the median n for galaxies 
with D/T > 0.5 being around 2.8±0.2 compared to 3.2±0.05 
for galaxies which are purely ellipticals ( D/T= < 0.1) The 
older galaxies also have slightly higher n values, and the de- 
pendence on n on Kirc is basically flat. 

We note that Hopkins et al. (2008) have recently shown 
that fitting a single component Sersic profile can mischar- 
acterize the profile shape of post starburst ellipticals. As 
we have seen with the profiles plotted in Section 3.2, kine- 
matically selected early types can have low Sersic indices. 
Hopkins et al. used two component Sersic fits to examine 
the inner and outer profiles of simulated merger remnants. 
A median Sersic index for the outer profile found by Hop- 
kins et al. is n ~ 2 — 3 (e.g. see their Figure 15). Although 
our mass and spatial resolution is not high enough to permit 
two component fits here, this issue should be borne in mind 
when using n to separate galaxies into morphological classes 
(e.g., see Section 4.3) 

The specific star formation rate, in the top row of Fig- 
ure [S] is pretty flat in the median with respect to the other 
properties. The strongest apparent trend is with galaxy age 
where it can be seen that both the median sSFR and the 
10%-90% envelope do drop with increasing age. The en- 
velopes do tend to widen also at low Vctt and Kirc, indicating 
that galaxies are more extreme with respect to sSFR when 
they are small. 

Looking at Figure [6] overall, we can see that there are 
many instances of the expected trends between properties 
which we know from observations (e.g., Blanton et al. 2003). 
For example galaxies are usually classifled into "early" and 
"late" types where one can broadly say that early types are 
larger, have older stellar populations and larger bulges than 
late types. Of the 15 panels in Figure |6] which carry indepen- 
dent information, 10 of them fit into this simple picture (for 
example, galaxies seen kinematically to have larger disks do 
have smaller Sersic indices), 1 runs counter to it (we find 
that older galaxies are smaller, as measured by reff) and 4 
are neutral (sSFR exhibits little median variation with other 
properties). It is beyond the scope of this paper to examine 
any of these relationships in detail, but it is apparent that 
simulations are nearing the point where in future work it 
will be profitable to do so and make more direct compar- 
isons with observational data. 

3.5 Angular momentum 

Simulations of disk galaxy formation have historically had 
problems with excessive angular momentum loss, leading 
to stellar components with bulges too large compared to 
observations (Mayer et al. 2008 for a review). These in- 
clude for example the zoomed simulations of Sommer-Larsen 
et al. (2003), Abadi et al. (2003ab), Robertson et al. (2004). 
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Figure 7. Rotation speed vs specific angular momentum j at 
redshift 2 = 1 for late type galaxies in the BHCosmo simulation 
(points). Only galaxies with a D/T ratio > 0.2 are shown, the 
same threshold used to define late types in the rest of the paper. 
Prom top to bottom are results for the stellar component, the gas 
and the dark matter. The bottom plot shows an estimate of the 
specific angular momentum that is not computed directly from 
the angular momentum of the particles but instead by multiplying 
the disk scale radius by twice the circular velocity. This is the 
estimator used for observational data (e.g., Courteau 1997) and is 
valid for an exponential disk. The solid line in each panel is taken 
from Abadi et al. (2003a), and describes well the median relation 
seen in the observational (z = 0) data compilation of Navarro 
(1998). In order to compare directly to Abadi et al. (2003a) we 
have rescaled using h = 0.65 when computing the values of j. 



Robertson et al. (2004) and Abadi et al. (2003ab) pointed 
out that exponential galactic disks could be formed which 
were of the correct size, but which were deficient in specific 
angular momentum. Governato et al. (2007) and Ceverino 
& Klypin (2007) have most recently been able to form disks 
which were closer to being realistic, through a combination 
of high resolution and modelling of star formation and feed- 



back. All of these simulations were zoomed runs, however, 
and it is not clear to what extent these results are dependent 
on the histories of the individual halos selected for resimu- 
lation. 

A clear way to describe the extent of any potential angu- 
lar momentum problem is to plot the specific angular mo- 
mentum content of galaxies as a function of their circular 
velocities. This has been done by for example Navarro & 
Steinmetz (2000) and Abadi et al. (2003a), comparing the 
simulation results against the observational compilation of 
Navarro (1998). Observations of the specific stellar angular 
momentum content of disk galaxies were made using an es- 
timator which combines the observed circular velocity and 
disk scale length (Courteau 1997), so that jobs = 2i?cft Vciic- 
For the simulations, the circular velocity is usually computed 
using Equation 2, but the angular momentum is summed di- 
rectly from the star particles. 

In Figure [7] we show results from the BHCosmo simu- 
lation at 2 = 1. It is important to note that the simulated 
galaxies we are plotting are from z = 1 whereas the observa- 
tions (Mathewson et al. 1992, Courteau 1997) were made at 
z ~ 0. Based on the lack of evolution in angular momentum 
vs Vcirc seen in the D4 run (which was run to z = 0) we 
expect this to make little difference to the comparison. The 
straight line in Figure [7] is taken from Abadi et al. (2003a). 
It describes the mean of the observational data well. For 
the simulation points, because we are comparing to obser- 
vational data for disk galaxies, we apply a threshhold on the 
D/T ratio before plotting points in Figure [T] This threshold 
is at a D/T ratio of 0.2, which selects 40% of the galax- 
ies to be disks at this redshift, z = 1. Although this cut is 
essentially arbitrary (we have a smooth continuum of D/T 
ratios for galaxies) , we will see later (Section 4.3) that it 
enables us to reproduce very approximately the fraction of 
morphologically classfied galaxies of early and late types as 
a function of density. For the points in Figure [7] we measure 
j, the specific angular momentum by summing the total an- 
gular momentum of the particles in question (we deal with 
stars, gas and dark matter separately in the different panels) 
and dividing by their summed mass. We use only particles 
within 2.2 times Res, in order to be consistent with ob- 
servations, and include all particles interior to that radius, 
including those in the bulge. 

We can see in the top panel of Figure Qthat the median 
jstar is a factor of ~ 8 lower than the observational value for 
disk galaxies. As we have seen in the distribution of D/T ra- 
tios, simulated galactic stars have too little orbital angular 
momentum. This result is somewhat better than the lower 
resolution runs of Navarro & Steinmetz (2000), who were a 
factor of ~ 30 below the observations, but the angular mo- 
mentum deficit is still large. Robertson et al. (2004), with 
somewhat better resolution in their resimulated galaxies (a 
dark matter particle mass of 3 x lO^h'^M©, 4 times smaller 
than Abadi et al. 2003a) formed a spiral with j'star approx- 
imately a factor of 2 below the observational trend. This 
galaxy, with a quiet late merger history is similar in its jatar 
to many of our galaxies. Our results are therefore likely to 
be similar to theirs, except that we have access to a whole 
population of simulated galaxies. 

The median gas specific angular momentum in Figure 
[7] is closer to the observational line than the stars. This was 
also seen by e.g., Abadi et al. (2003a) and Robertson et al. 
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Figure 8. The stellar mass TuUy-Fisher relation for late type 
galaxies (selected using D/T > 0.2) in the BHCosmo simulation 
(at redshift 2 = 1 and 2 = 3) and the D4 simulation (z = 0). 
The X-axis is the statistic So. 5 from Kassin et al. (2007), which 
is a kinematic estimator of support from disordered and ordered 
motion (see equation 3). Sq.s from the simulations was computed 
either using gas particle velocities (left panels) or star particle 
velocities (right panels). The solid line shows a best fit power law 
relation to the simulation data in each panel. The dashed line 
is the best fit power relation from Kassin et al. (2007) to their 
observational data. This observational fit is the average relation 
measured from 4 bins varying in redshift from z = 0.1 — 1.2. 
Kassin et al. (2007) find no statistically significant variation with 
redshift. 



(2004), who find that as gas cools and stars form, angular 
momentum is lost to the galactic halo. For the dark matter, 
the ji'dm values are closer to those seen in the stars, (as also 
seen by Abadi et al. 2003a) , although the scatter in jdm for 
a given Vdrc is larger. 

Although the angular momentum content of the disks is 
low, we find that the size of the disks is actually reasonable 
compared to observations. To show this, we plot in the bot- 
tom panel of Figure [7] the observational estimator for the 
specific angular momentum of stars, jobs- In this case, we 
find that the median is compatible with the observational 
line. Therefore, as seen by Roberton et al. (2005) in a single 
galaxy, the disks are of reasonable scale, but they have a 
bulge component which is too large, so that too much of the 
galaxy is supported by non-circular motions. 

By looking at the angular momentum content of ~ 500 
disk galaxies, we see that the problem seen in resimulations 



of a small number of galaxies persists. The addition of black 
hole accretion and feedback to our runs compared to those 
mentioned above does not noticeably impact the excess pop- 
ulation of non-rotationally supported stars. This is compat- 
ible with the recent work of Okamoto et al. (2007) who res- 
imulate disk galaxies with and without black hole feedback. 

3.6 Stellar mass Tully-Fisher relation 

The Tully-Fisher (1977) relation (hereafter TFR) between 
the luminosities of galaxies and their rotation velocities has 
long been used as a test of galaxy formation models (e.g., 
Steinmetz & Navarro 1999). For example, the amplitude and 
zero point of the TFR seen in observations at 2 ~ is 
still difficult to reproduce in high resolution simulations of 
disk galaxies (Portinari & Sommer-Larsen 2007). The stel- 
lar mass-rotation velocity version of the TFR is easier to 
compute theoretically. It has recently been measured out 
to z — 1.2 using observational data from the AEGIS and 
DEEP2 surveys (Kassin et al. 2007). In the usual M-Vrot 
version of the TFR, where M is the galaxy stellar mass and 
Vrot there is a strong dependence of scatter in the relation 
on galaxy morphology (disturbed, close pairs etc.). Kassin 
et al. (2007) however find a modification of the TFR which 
works well with minimal morphological pruning and holds 
over the entire redshift range z — 0.1 — 1.2. 

The relation computed by Kassin et al. is between M 
and So. 5, where the latter is a kinematic estimator which 
combines dynamical support from ordered motion with that 
from disordered motion (see Weiner et al. 2006). The quan- 
tity So.5 is computed using 



C.2 
■50.5 



0.5Vrot^ + CTg, 



(3) 



where Vrot is the rotation velocity of galaxy on the flat part 
of the rotation curve, ag is a gas velocity dispersion, both 
being measured from a two-dimensional fit to emission line 
in the galaxy spectra. Kassin et al. find the observed M — So.s 
relation to be fairly tight, with an rms intrinsic scatter in 
So. 5 of 0.10 dex. The slope and amplitude of a log- log fit 
to the points, logj^Q So.r, = a + blogjQ AZ/lO^^h^^Mo is a = 
0.34±0.05 and b = 1.89±0.03 averaged over 4 bins in redshift 
ranging from z = 0.1 — 1.2. The fits are not significantly 
different for all redshift bins. 

We compute M and So.s values for disk galaxies in our 
simulations. Disk galaxies are again taken to be those for 
which D/T > 0.2. The M values are straightforward: we use 
the total mass in star particles in each galaxy's subhalo. For 
So. 5 we have a choice, whether to use a quantity measured 
from particle velocities in the simulation, or whether to use 
the mass within a certain radius to compute a value of Vrot, 
assuming rotational support. In order to include the effect 
of non-rotational motions, we do the former. For Vrot, we 
compute the mean rotational velocity of particles in a ra- 
dial bin of width rcfi centered on a radius 2.2rcff from the 
center of each galaxy, and we compute ag from the scatter 
of particle velocities about V-ot- Using these quantities, we 
estimate So. 5 for the simulated galaxies. For each galaxy, we 
compute two versions of So. 5, one using gas particle infor- 
mation, which more closely corresponds to the observational 
information, and another using the star particles. 

The results are shown in Figure [S] where we show So. 5 
computed from both gas and stars for 3 different redshifts. 
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For 2 = 1 and 2; = 3, we use the BH Cosmo simulation, 
whereas at 2: = we use the D4 run. We also show the 
power law relation which is the average one for the data 
of Kassin et al. (2007) over 2 = 0.1 — 1.2. For each panel, 
we also compute the best fit power law relation from our 
simulated data and show it as a solid line. In computing the 
fit, we follow Kassin et al. and assign the intrinsic scatter (in 
this case we assume all the scatter is intrinsic) to the 5*0.5 
coordinate. 

We notice from Figure |8]that the 5*0.5 values computed 
from the gas have a larger scatter about the mean relation 
than the stars. The rms scatter in So. 5 ranges from a max- 
imum of 0.16 dex for the gas at 2 = to a minimum of 
0.08 for the stars at 2 = 1. This is not to dissimilar from 
the observational data, which has So. 5 instrinsic scatter ran- 
gin from 0.08 to 0.12 dex over 4 redshift bins. As with the 
observations, we find no systematic change in rms scatter 
with redshift. It can also be seen that fits to the stellar 
simulated So. 5 values are consistent with the observed fit. 
The largest discrepancy is at 2 = 3, a higher redshift than 
yet reached with the observational data, for which the slope 
is 6 = 0.24 ± 0.02, less than 2a different. The So. 5 values 
computed from the gas particles do appear to have a sys- 
tematically slightly steeper slope than the observations (we 
find b = 0.22 ±0.03 averaged over all 3 redshifts). but this is 
less significant than the ~ 4(t discrepancy in the intercept, 
a which is seen at 2 = 0. 

Overall, we therefore find that the stellar mass TFR 
relation (between M and So. 5) measured in the simulations 
does not change in slope with redshift (between 2 = and 
2 = 3) and is consistent with the slope seen in the observa- 
tional data of Kassin et al. (2007). The zero point measured 
from the gas properties in the simulation is somewhat noisy, 
but not systematically different when averaged over all red- 
shifts. 

3.7 Orientation 

The relationship between a galaxy's morphology and angu- 
lar momentum can give important clues to the processes 
that occurred during its formation. A way to measure this 
quantitively is to measure the alignment between the first- 
ranked principal axis of intertia of a galaxy and its angular 
momentum axis. We have carried this out for all galaxies in 
the BHCosmo simulation, again at 2 = 1, including the low 
D/T ratio galaxies. 

We compute the inertia tensor for each galaxy and di- 
agonalize it, to find the principal axis with the largest mo- 
ment of inertia. We do this separately for each component 
of the galaxy, using every particle in the subgroup. From 
the scalar product of this axis with the angular momentum 
axis, we find an alignment angle, 6. The pdf of 9 values is 
shown in Figure [Q] We find that the gas component of galax- 
ies is most strongly aligned, with a pronounced peak close 
to zero degrees, likely to be coming from the gas disks seen 
in many galaxies (for example, even the early type galaxy 
in the top panel of Figure[2]had a pronounced gas disk. The 
median value for the gas is 6 = 25.7°, significantly more 
aligned than the stars (9 = 39.0°). The stellar alignment 
angles do have a low 6 peak, however, indicating that rota- 
tionally supported stellar disks are still common. By looking 
at the orientations of galaxy shapes in Figure [2] (where the 
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Figure 9. Histogram of the angle (in degrees) between tlie prin- 
cipal axis measured from the inertia tensor and the angular mo- 
mentum axis, for galaxies in the BHCosmo simulation at redshift 
2 = 1. We show results for the stars, gas and dark matter distri- 
butions, but the angle in each case is measured from the angular 
momentum axis computed from the stars. The solid points show 
the median angles (25.7° for gas, 39.0° for stars, and 48.3° for 
dark matter). 

projections are defined by the A.M. axis) we can see that 
the misaligned galaxies are indeed likely to be mostly early 
types, as we would expect. Finally, in Figure |§] we can see 
that the dark matter halos are also noticeably flattened per- 
pendicular to the rotation axis (the same flattening direc- 
tion as the gas and stars), with a median alignment angle 
oi 9 = 48.3°. The expectation value for the median of a 
random distribution is 9 = 60.4°, with a Poisson error of 
1.8°. As the dark matter is dissipationless, this flattening 
could be partly due to adiabatic contraction in response to 
the rotationally flatted gas and stellar component. A strong 
effect was also seen by Bailin & Steinmetz (2005) in purely 
dissipationless simulations. 

Looking at the alignments of the angular momentum 
vectors of the different components themselves can also give 
us some clues to how the rotational support arises. In fig- 
ure [TD] we plot the distribution of alignment angles (f> of the 
star and gas angular momentum axes with the dark mat- 
ter angular momentum axis for each galaxy. We find that 
there is a strong alignment of the stellar A.M. with the dark 
matter A.M. (median angle = 43.5°), and a somewhat 
weaker alignment of the gas with the dark matter (median 
(f> — 69.6°). This is quite suprising given the alignment of 
the principal axes of the gas seen in Figure (9] We have also 
computed (not plotted) the alignment of gas and star an- 
gular momentum vectors. We find a median angle in this 
case of (/) = 61.6°, closer to the alignment of the gas and 
dark matter. We can therefore see that the stellar angular 
momenta appear to have become correlated with that of the 
dark matter after the stars have formed from the gas. This 
coupling through coUisionless dynamics could perhaps be a 
numerical effect, due to two body interactions. To provide 
additional information, we have investigated our lower res- 
olution simulation, the D4 and found that the rank order 
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Figure 10. The alignment angle, (in degrees) between the 
star and gas angular momentum axes and the dark matter angu- 
lar momentum axis, for galaxies in the BHCosmo simulation at 
redshift z = 1. The solid points are the medians of the distribu- 
tions (43.5° for stars and 69.6° for gas). 

of the median alignment angles are the same as for the BH- 
Cosmo run, i.e. star-dark, star-gas, gas-dark, in order of best 
to worst alignment. The median angles of the latter two are 
~ 10° worse, but for star-dark, the angle is ~ 20° worse, 
suggesting that the stellar angular momentum is affected by 
the most by lack of resolution. 

We note that in the simulations of Van den Bosch et al. 
(2002), a much stronger alignment was seen between the 
angular momentum of the (adiabatic) gas in halos and the 
dark matter (median angle cj> — 27.1°) than here, so that 
without dissipation, the gas is much more strongly coupled 
in its rotation to the dark matter. 



3.8 Age of disk and bulge stars 

We have seen in Figure [6] that galaxy ages do have some de- 
pendence on other properties, including size and D/T ratio. 
We have decomposed the galaxies kinematically into bulges 
and disks (Section 3.1), and so it is of interest to see whether 
the ages of the disk and bulge stars fit into this general pic- 
ture. In Figure 1111 we have plotted the mean redshift of 
formation of stars in the disks and bulge components of the 
galaxies in the BHcosmo run at z = 1. We show the galax- 
ies with D/T ratios greater than 0.2 (our fiducial threshold 
between kinematically defined early and late types) as red 
points. 

We can see that the cloud of points extends down to 
low redshifts, indicating that the disk component of galax- 
ies does indeed form significantly later than the bulge. For 
example, there are essentially no bulges with a mean star 
formation redshift below z = 3, but disks continue to form 
all the way down to z — 1, the lowest redshift possible. The 
disks below « ~ 2 are all disk components of bulge domi- 
nated galaxies, so that they only represent a small fraction 
of the galaxy. Interestingly, we also see that the disks with 




disk formation z 

Figure 11. The mean age of star particles defined kinematically 
to be in a disk component, versus the age of those in the bulge, 
for galaxies in the BHCosmo simulation at redshift z = 1. We 
show results for late type galaxies (D/T > 0.2) as red points and 
early types as black points. 

the highest formation ages z ^ 4.5 also belong to early-type 
galaxies. 



4 RELATION BETWEEN ENVIRONMENT 
AND GALAXY PROPERTIES 

In the previous section we have seen that individual galaxy 
properties in the simulations are by and large what we would 
expect given a mixture of early and late types. We will 
now investigate how the immediate environment has affected 
these properties. 

4.1 Measures of environment 

We use four different definitions of the local environment of 
a galaxy, in order to span a range of parameters from those 
close to what can be measured observationally to measures 
related directly to the dark matter distribution. The first 
measure is the density computed from the distance to the 
tenth nearest galaxy. For this we use positions for galax- 
ies above our 5000 particle threshold, and compute a local 
density for each galaxy, dividing the density of the sphere 
enclosing 10 neighbours by the mean density of galaxies 
in the simulation. We label this quantity on the plots as 
p/ (p). The density computed from the tenth nearest galaxy 
is often used in observational measurements, such as in the 
morphology-density relation of Dressier (1980), and subse- 
quently by Dressier et al. (1997), Smith et al. (2005), Carpak 
et al. (2007) and others. In these cases, the measurement is 
made in 2 dimensions, from the projected positions of galax- 
ies, but the two quantities will be closely related. We explore 
this measure more closely in the next section. 

A related measure is pis, which is the density com- 
puted from a sphere centered on the galaxy which contains 
10^'^h~^MQ of dark matter. In practice, this will be very 
similar to the first measure, based on galaxies, but will be 
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less noisy. It will also be more directly related to the gravi- 
tational influence of nearby structures, being based on mass 
density and not number density. We note that a very sim- 
ilar measure of environment, computed from the radius of 
a sphere which contains lO^^h~^M0 of all types of matter 
was used by Colberg & Di Matteo (2008) in their study of 
the environmental dependence of black hole properties in 
the same simulation (the BHCosmo run) which we use here. 

The galaxies are members of friends-of-friends groups 
of particles, and so we can also categorise their environment 
using these grops. We use the mass of the FOF group as 
one means of quantifying their environment, and the other 
the distance of the galaxy from the FOF group center of 
mass, in units of the virial radius, 7-200 • We compute r2oo by 
measuring the radius of a sphere which contains 200 times 
the mean density. 

4.2 Galaxy properties and environment at z=l 

In Figure [T2I we use the same galaxy parameters which were 
shown in Figure [6] but this time instead of plotting them 
against each other, we plot them against our four measures 
of local environment. We again show the situation at z = 
1, for the BHcosmo simulation. Before proceeding with our 
analysis we note that ai z — 1 there are no rich clusters 
in the relatively small simulation volume, (the largest FOF 
group has a mass of 6 x lO^'^h^^M©). Below we will see that 
this will affect our ability to compare to observational data 
for the high density end of environments at this redshift. 
In interpreting Figure [T^] we should therefore bear in mind 
that the dense environments in the plot are groups rather 
than clusters. In order to look at galaxies in clusters at high 
resolution, a resimulation approach should be used (see for 
example Saro et al. 2006.) 

Looking at the figure, we can immediately see that the 
results are a mixture of relationships which we would have 
expected and some which are surprising. The values of Vcirc 
for galaxies increase steadily with density, both p/ (p) and 
pi3, the median rising by approximately a factor of 2 over 
3 orders of magnitude in density. For Vdvc vs Mgrp, we see 
that there is an upper envelope of large galaxies in each 
group, but also that when Mgrp reaches 2 x lO^'^h~^M0, 
there are enough satellite galaxies to pull the median Vdic 
value down. In terms of r/r2oo, we can see that the galaxies 
with the very highest Vcirc values are indeed close to the 
center of mass of the FOF groups. There is a large cloud of 
points, between 0.1 and ~ 1.0 virial radii from the center, 
indicating where most of the galaxies are found. The number 
density of points decreases rapidly beyond r/r2oo, although 
a few galaxies are found at distances of up to 10 virial radii. 
The median of Vdic values does not vary with r/r2oo. 

The variation of r^g with density is similar to Vdic, 
except that the effect of satellite galaxies appears to be more 
pronounced. Also, even the median r^s value changes with 
r/r2oo. 

Moving on to the D/T ratio in the third row from bot- 
tom of Figure 1121 we find more complex dependencies on 
density. There is a moderate rise in the median D/T ratio, 
from 0.13 to 0.2 as p/(p) > changes by 3 orders of mag- 
nitude. It therefore seems as though there are many small 
bulge-heavy galaxies at low densities. This behaviour is also 
repeated but more weakly (only significant at the ~ la level) 



when looking at pia. If we look at the relationship between 
D/T and r/r 200 we can see that the median D/T rises to- 
wards the edges of groups, significant at the 3 — 4(t level, 
and the 90th percentile line increases dramatically. There 
are therefore not many disk dominated galaxies in the cen- 
ters of groups. This seems to be opposite to the trend seen 
in the variation of D/T with density. However, it is possible 
that the D/T values do go down with very high densities, but 
that this is not easy to see with the relatively noisy density 
estimators that we have (the 90th percentile of D/T does ap- 
pear to decrease for large P13.) When comparing r/r2oo and 
p values we should bear in mind that when large galaxies 
are present, the group center of mass better correlates with 
the centers of large galaxies, and hence the density. This can 
help explain some of the apparent differences in behaviours 
with r/r2oo and p. 

The variation of Sersic index, n with density is basically 
consistent with flat from low to moderate density, followed 
by a decline in the median value, indicating fewer centrally 
concentrated early types at high density. This is somewhat 
consistent with the trend seen in D/T ratios, although in 
that case more variations were seen at low density than high. 
For n, the errors are large, but it does seem from the 90th 
and 10th percentiles that late types preferentially prefer the 
outskirts of groups as well. 

The mean stellar age of galaxies does not appear to vary 
much with density at this redshift. There is perhaps even a 
small hint that galaxies in low density regions are slightly 
older, something which is puzzling but consistent with the 
slight preponderance of bulge-dominated systems there. 

The specific star formation rate of galaxies in our sim- 
ulations rises slightly with density, but also on the outskirts 
of groups. While the trend seen in observations at 2: = 1 by 
Cooper et al. (2008) has the same sign as that seen at low 
redshifts (e.g., Gomez et al. 2003), Cooper et al. find that 
the total SFR has an inverted relation with density. The be- 
haviour and comparison with observations, particularily of 
the total SFR is investigated fully in this BHCosmo simula- 
tion by Colberg & Di Matteo (2008) . For now, we note that 
the specific SFR trend we plot is roughly consistent with 
what has been seen with the D/T ratio and Sersic index, 
i.e., a reversal, or at least flattening of the trends seen in 
galaxy morphology with environment at redshift z — 0. 

Some of the trends of galaxy properties with density 
do have some dependence on galaxy mass/particle number. 
For example when we limit ourselves to plotting galaxies 
containing 25000 particles instead of 5000 (although we still 
use all galaxies to compute p/ (p)), we flnd that the me- 
dian D/T value for the highest p/ (p) bin drops by 2 cr to 
D/T=0.1. The previous bins are essentially unchanged. It is 
therefore possible that the efficiency of destruction of spirals 
in the highest density regions depends on particle number. 
In this test, the qualitative behaviour of SFR is unchanged, 
but the highest density bin does have slightly older galaxies. 

Our overall conclusion from Figure [12] is that at 2; = 1 
the relationships between galaxy properties and density in 
our simulations are complex. The preponderance of late 
types in low density regions seen at z = seems to be 
reversed, with higher star formation rates in high densi- 
ties. At the very highest densities, there are hints that the 
disk fraction may decrease, but this is dependent on reso- 
lution/particle number, with the main uncertainty that our 
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Figure 12. Galaxy properties versus environment at 2 = 1 in the BHCosmo simulation. The measures of environment are, from left 
to right, galaxy density (measured from the distance to the 10th nearest neighbour galaxy) in units of the mean, mass of the host 
friends-of- friends group, distance from the center of mass of the friends-of-friends group, in units of the virial radius, rgoOi ^■nd dark 
matter density computed from the radius which encloses lO'^^h^^MQ of dark matter. In each panel, results for individual galaxies are 
shown as points, and the median is shown as a thick solid line with Poisson error bars. The 90th and 10th percentiles of the distribution 
in each bin are shown as thin solid lines. 



simulation box, at 33.75 /i~^Mpc is too small to contain any 
clusters at z = 1. Cooper at al. (2007) have investigated the 
observed relationship between galaxy color and local den- 
sity from redshifts z — 0.4 to 2 = 1.35, finding that the 
trend of red sequence galaxies to favour overdense regions 
evolves strongly over that redshift range. The locations of 
red galaxies z ^ 1 are however still seen to correlate with 
high densities. Population synthesis models will be needed 
to compare colors of our simulated galaxies with this infor- 
mation. For now, we restrict ourselves to the morphological 



information we have computed, although as we shall see in 
the next section, the morphology-density relation in the sim- 
ulation also appears to evolve somewhat more weakly than 
in the observations. 



4.3 Morphology-density relation 

The tendency of dense (galaxy cluster) environments in the 
low redshift Universe to contain a higher fraction of early 
types than the field has long been known (Dressier 1980). 
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Figure 13. The morphology-density relation at 3 redshifts. We 
show as different lines results for galaxies types selected on the 
basis of kinematic D/T ratio, Sersic index, and mean stellar age 
(see text). The z— axis, S is the projected galaxy density. Panels 
for z = 1 and ^ = 3 are for galaxies from the BHcosmo simula- 
tion, and z = galaxies are taken from the D4 run. As shaded 
histograms in the top 2 panels we show the observational results 
of Dressier (1980,1997) at z ~ and Smith et al. (2005) at 2 = 1. 



Recently it has become possible to extend this work to 
higher redshifts, z = 0.5 — 1 and beyond. Here we compare 
simulation results for this morphology-density relation to 
the observational data of Smith et al. (2005), which, along 
with data from z = 0.75 — 1.27 includes a compilation at 
lower and intermediate redshifts from Dressier (1980, 1997), 
and Treu et al. (2003). 

The galaxies in the Smith et al. sample and compilation 
have limiting V magnitude of -20.2 at 2 = (-21.2 at z=l). 
Their space density (using the luminosity function of Brown 



et al. 2001) is 5.12 x 10"^ ( ft"^Mpc)"^ at z=0. The observa- 
tional fraction of early-type galaxies is plotted as a function 
of projected galaxy density at 2 different redshifts in the 
top 2 panels of Figure 1131 The projected galaxy density is 
the density enclosed within a rectangular region extending 
to the 10th nearest galaxy. In order to include redshift in- 
formation in the z — 1 data. Smith et al. include galaxies 
in the density estimate in slices of redshift width Az = 0.1. 
They also subtract a background galaxy density from the 
density estimate and also multiply by a density correction 
factor of 2. After this procedure, their results are that the 
correlation between early type fraction and density is still 
present at z = 1, but significantly weaker than at 2 = (and 
also 2 — 0.5, which lies in between, and is not plotted here). 

In order to make predictions from the simulations, we 
also compute a projected galaxy density measured from the 
area delimited by the 10th nearest neighbor. It is not pos- 
sible to replicate entirely the relatively involved procedure 
employed by Smith et al. , but we make sure that we are cen- 
tered on the same mean projected galaxy density, and make 
sure that our results are robust by also trying a completely 
different density estimator, as explained below. 

We project the whole simulation box, using information 
from the 3 different orthogonal projections in our data sam- 
ple. This results in a lower mean projected galaxy density 
than Smith et al. , who at z = 1 have a mean area density 
of galaxies of 5.25 per Mpc^. Our galaxy density is 2.0 per 
Mpc^, so that we correct our projected simulation density 
values by a constant factor to bring them into agreement. 
This is under the assumption that the overdensity of galax- 
ies relative to the mean is the physically relevant quantity. 
We have also tried computing the morphology density re- 
lation using the 3 dimensional galaxy overdensity from the 
simulations, for comparison. We find that all the qualitative 
conclusions we draw from the comparison with projected 
observational data also hold, giving us confidence in the ro- 
bustness of the results. 

In order to find the early type fraction, we compute 
three difi^erent values using three different criteria, kinematic 
D/T ratio, Sersic index and age. For the kinematic selection, 
we use a cutoff of D/T=0.2 between early and late types, 
and for the Sersic index, a threshold of n = 3.0, For the age 
criterion, we used a cutoff mean stellar formation redshift 
between early and late types which is different for the 3 red- 
shift bins (2 = 0.0, 1.0 and 3.0) we use. In order to set the 
threshold, we choose values which make the overall fraction 
of early and late types the same as for our D/T ratio cri- 
terion. This cutoff stellar formation redshift is 2 = 2.6, 3.8 
and 4.8 for redshift bins 2 = 0, 1.0 and 2 — 3.0 respectively. 

Our results for the morphology-density relation in the 
simulations are shown in Figure [TH] alongside the observa- 
tional results at 2 = and 2 = 1. We also show the simula- 
tion predictions for redshift 2 = 3. We can see that at 2 = 0, 
the early type fraction increases in high density regions. Al- 
though there are no rich clusters in our relatively small vol- 
ume, the simulation curves with early types selected kinet- 
matically by and age track reasonably well the observational 
results up to moderate galaxy density. The photometric cri- 
terion however diverges at the high end of the density range, 
with the number of galaxies with n > 3 falling dramatically 
at projected densities > 20Mpc~^. This behaviour was also 
seen in Figure [T2l 
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At redshift z = 1, the morphology-density relation in 
the simulation flattens out, with no increase seen in the frac- 
tion of early types at high densities in the middle panel of 
Figure 1131 This change between z = and 2 = 1 is some- 
what more extreme than that which is seen observationally, 
although the sign of the change is the same. Again the se- 
lection by Sersic index gives the most discrepant results. At 
z — 3, the simulations predict a complete reversal of the 
morphology density relation, with the fraction of galaxies 
with younger stellar populations being highest in the dens- 
est regions. It is to be noted that because the units of the 
cc-axis of Figure [T3] are in proper Mpc, the galaxy density of 
the different redshift samples is affected by the expansion of 
the universe. 

Because the simulation volume is relatively small, we 
are restricted to computing the morphology-density relation 
for galaxies with a relatively small range of stellar masses 
(there are few extremely massive galaxies for example, par- 
ticularly at 2 = 1 and above). Recently, it has been found 
(F. van den Bosch, private communication, Weinmann et al. 
2008) that for (satellite) galaxies in a fixed stellar mass bin, 
there is no significant morphology dependence on host halo 
mass, i.e. observationally the morphology-density relation 
seems to be driven by the stellar mass-density relation. The 
results we find in the simulation for the morphology-density 
relation could be therefore be relatively fiat at least partly 
because of our small range of galaxy masses. We note that 
we have accounted for differences between central and satel- 
lite galaxies in our study, whereas Weinmann et al. find that 
their result holds for satellite galaxies only. Larger simula- 
tions will be needed to investigate conclusively whether the 
morphology-density relation in simulations is fiat for bins of 
galaxy stellar mass. 

We have seen that the morphology-density relation flat- 
tens out or inverts at higher redshifts. However, what is not 
easily seen in Figure [13] is that the overall fraction of early 
type galaxies change with redshift as well. We have applied 
the same thresholds to separate early and late types (D/T 
> 0.2 and n < 3.0 for late types) in outputs of the BHCosmo 
run at redshifts up to 2 = 9. The results are shown in Fig- 
ure [141 where we can see that the two means of separating 
galaxy types do agree at 2 = 4 and below. At high redshifts, 
although the errors are large, the correspondence between 
type chosen through D/T and n is not good, with galaxies 
appearing to become steeper of profile but at the same time 
having smaller kinematically defined bulges. Based on the 
resolution tests (Figure [5}, D/T is likely to be a more ro- 
bust discriminator between galaxy types. This conclusion is 
understandable if we take into account the fact that single 
component Sersic profiles are often not a good character- 
ization of the profile shape of post-starburst ellipticals, as 
found by Hopkins et al. (2008) and discussed further in Sec- 
tion 3.4. 

At lower redshifts, however, the number of early types 
does increase significantly. As described by D/T ratio, the 
early type fraction increases from ~ 0.35 at 2 ~ 9 to ~ 0.6 at 
2 = 1. The early types chosen by n also increase below 2 = 4 
Following galaxy properties in detail with redshift is beyond 
the scope of this paper, so we will confine ourselves to noting 
that the general decrease in the fraction of disk galaxies as 
the universe evolves seen here is the same qualitative trend 
seen observationally (e.g.. Postman et al. 2005). 




Figure 14. The early-type fraction of simulated galaxies as a 
function of redshift. Galaxies were classified either on the basis of 
kinematic D/T ratio (> 0.2 for late type) or Sersic index, (n < 3.0 
for late type), as in Figure [l3]l. The filled points (i.e. z ^ 1) are 
for the BHCosmo simulation and the open points for the D4 run. 



5 LARGE-SCALE STRUCTURE AND 
CLUSTERING 

The galaxies in our simulations form at the same time as 
the cosmic web of large-scale structure. We have seen from 
the previous section how the local environment of galaxies 
affects their properties. We will now turn to larger scales, 
examining the morphology of this structure as traced by dif- 
ferent types of galaxies and also their clustering, measured 
using the autocorrelations function. As with the previous 
section, the relatively small simulation volume means that 
we will not have instances of truly rare events, rich clusters 
of galaxies, and also the box size will have some effect on the 
overall amplitude of clustering. For this reason, we will con- 
centrate on the differences between dark matter and galaxy 
clustering, their relative bias, for different galaxy types. This 
should be less sensitive to the small box size than the raw 
amplitude of the autocorrelation function (following the ap- 
proach of e.g., Katz et al. 1999). 



5.1 Visual impression 

In Figure [TS] we show 10 /i~^Mpc thick slices through the 
galaxy and dark matter density fields, at 3 different red- 
shifts. We have used our usual D/T theshold of 0.2 to clas- 
sify galaxies into early and late types, which are plotted 
using different symbols. The same 5000 particle lower limit 
for subhalos is again used at all redshifts, except for the left- 
most panel, which is taken from the D4 run. In this case, we 
use the equivalent particle limit, scaled from the BHCosmo 
run using the ratio of mass resolutions. 

Looking at the three different redshift panels, the over- 
all framework of structures looks quite similar at different 
times. In the dark matter plots, the dense regions at 2; = 3, 
are quite web-like, made up of many intersecting filaments, 
whereas at 2 = they have collapsed into more compact 
structures. On larger scales, however, not much appears to 
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Figure 15. Large-scale structure in the galaxy and dark matter distribution. Top row: slices (10 h~^Mpc thick) through the simulated 
galaxy distributions at 3 redshifts. Late type galaxies (selected so that D/T > 0.2) are shown as spiral symbols and early types as points. 
The leftmost panel, 2 = is for the D4 simulation, and the two right panels (z = 1, middle and 2: = 3, right) are from the BHCosmo run. 
Bottom row: a linear greyscale plot of the dark matter distribution in the same regions of the same simulations shown in the top row. 



have changed, particularly when looking at the galaxy dis- 
tribution. 

Comparing the distribution of early and late types, we 
see some evidence for dense clumps of early types at z = 
(for example at the bottom of the panel). There may also 
be marginally more late types around the outskirts of these 
structures, although the situation does not look markedly 
different at higher redshift. At redshifts 1 and 3, however 
there do seem to be many early types which are relatively 
isolated, somewhat more than late types. More convention- 
ally, the protogroup of the bottom of the panel appears at 
redshift 3 to consist of more late types, and then at z = it 
collapses into a virialized group of early types. Apart from 
this, however the visual impression at z; = 1 and above is 
that late types are marginally more highly clustered in the 
simulation. We shall see below that this is borne out by 
quantitative measurement. 



5.2 The autocorrelation function 

We use the same delimitations between galaxy types as in 
Figure[T3]to investigate the autocorrelation function of early 
and late types. We look at three redshifts, 2: = 3, 2: = 1 and 
z = 0, using the D4 simulation for the latter. In addition, we 
separate the full sample into two sets, by mass, in order to 
see whether the mass dependence has a larger effect on clus- 



tering. In order to compare the results directly with those of 
galaxy types, we chose to make the higher mass subsample 
in each case have the same number of galaxies as those with 
D/T < 0.2. 

We compute ^(r) in real space, making use of the peri- 
odic boundary conditions of the simulation. The results are 
shown in Figure [161 Also given in each panel are values for 
the large-scale bias between galaxies and dark matter. For 
the early types it is computed using &e = ■\/ /£.dm. and 
analogously for the late types. We also compute the bias for 
the high mass and low mass subsamples. In each case we 
give the mean value for the 5 bins between 2 h~^Mpc and 
10 /i~^Mpc, and an error which comes from the standard 
deviation of the results for the 5 bins. 

On scales r 1 h~^Mpc we can see that in all pan- 
els we have approximately scale independent bias between 
galaxies and dark matter. All galaxy subsamples at z ^ 1 
have a higher amplitude of clustering on these scales than 
the dark matter, with the lowest amplitude corresponding 
to 6 = 1.08 ±0.03. At redshift z — 3, this bias can be as high 
as 6 = 2.48 ± 0.12, so that the galaxies follow the familiar 
"high peaks" biasing trend where similar mass galaxies at 
higher redshifts are associated with rarer, higher peaks and 
are biased more strongly (e.g., Bardeen et al. 1986). The val- 
ues of bias decrease until at 2: = 0, all galaxies are antibiased 
with respect to the dark matter. 
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Figure 16. The autocorrelation function of simulated galaxies as a function of type and redshift. The 3 columns show results for 2 = 0, 
2 = 1 and z = 3, from left to right, with the 2 = results taken from the D4 simulation and the rest the BHCosmo run. In each row, 
(a)-(c) we show galaxies classfied into early and late types by D/T ratio, Sersic index and mean stellar age. In row (d), we show results 
for galaxies split into two mass classes (see text). The galaxy correlation function results are shown as dotted lines for late types (and 
low mass, "small" galaxies in row (d)) and dashed lines for early types. In each panel we also plot the autocorrelation function of the 
dark matter as a thin solid line. The bias factor between the dark matter and the galaxies is given in each panel, computed from the 
average of the bias factor in the 5 bins with 10 h^'^Mpc > r > 2 h~^Mpc. The quoted 1 a errors are derived from the standard deviation 
of results for the 5 bins. 



The subsamples based on mass also act in the fashion 
expected from high peaks theory, with the more massive 
galaxies being more biased. However, if we look at the dif- 
ferent methods for classifying types, we can see that in many 
cases the late types are more highly biased than the early 
types. At z ^ 1, no matter how galaxies are classified, the 
late types are more highly biased. The largest difference is 
for the kinematic classification, at 2 = 1 for which high 
D/T galaxies are 25% more biased. Even for galaxies clas- 
sified based on mean stellar age we see this trend, with the 
youngest galaxies (assumed to be late types) being 7% more 
biased at this redshift. 

We therefore find that at high redshifts, the disk galax- 
ies have a higher correlation function (6 = 1.35 ± 0.03 at 
z = 1) on large scales than the high mass galaxy subsam- 



ple {b — 1.25 ± 0.05 at 2 = 1), even though they actually 
have a lower median circular velocity (iidrc = 136 kms~^ 
vs Ucirc = 177 km s"'^). This is quite interesting, implying 
that the processeses which determine whether a galaxy is 
an early or late type at these redshifts are even more closely 
coupled to clustering than galaxy mass. 

At redshift zero, the situation reverses itself for the most 
part, with early types chosen on the basis of D/T and mean 
stellar age being more clustered than the late types. The dif- 
ference in the clustering of the large and small mass subsets 
also becomes statistically insignificant. 

If we now turn to small scales r ^ 1 h~^Mpc, the 
clustering between galaxies and mass does become scale- 
dependent. Selecting galaxy type using n changes the am- 
plitude the most, with early types being a factor of ~ 20 
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times more clustered at a separation of r = 0.1 /i"^Mpc at 
2: ^ 1. This can be compared with the D/T plot, for which 
kinematic separation has little effect on ^. The difference in 
the slope of the stellar profile fitted for galaxies with different 
n therefore has effects on much larger scales than individual 
galaxies. Given that the effects of galaxy exclusion are seen 
in the large mass subsample ai z — 1 (^large drops to zero at 
radii presumably close to twice the radius of large galaxies) , 
one might worry that this is also happening in the n panel. 
However, looking at Figure [6] one can see that that galaxies 
with n > 3 essentially sample the entire Vdic range, so this 
is unlikely to be a problem. 

Observationally, several recent studies have begun to 
probe the details of clustering as a function of galaxy prop- 
erties at redshift z — 1. For example. Coil et al. (2006) 
have investigated luminosity dependence of clustering in the 
DEEP2 redshift survey, finding that more luminous galaxies 
are more strongly clustered than fainter ones, with the de- 
pendence of relative bias on luminosity being even stronger 
at z = 1 than z — 0. For a concordance cosmology, the bias 
seen by Coil et al. varies from b — 1.26 — f .54. The number 
density of galaxies for the observed sample with the low- 
est luminosity threshold is 1.3 x fO~^( h~^Mpc)~^ , within 
10% of the space density of our simulation sample of "large" 
galax;ies plotted in panel (d) of Figure[T6l The observed value 
of the bias from Coil et al. is b = 1.26 ± 0.04, compared to 
our simulation value of biarge = 1.25 ± 0.05. The observed 
clustering amplitude for a given space density therefore ap- 
pears to be reproduced by the simulations. We note that in 
this comparison, the bias measured by Coil et al. is the linear 
bias whereas we measure the bias with respect to the non- 
linear dark matter. At this relatively high redshift and on 
the scales we measure we do not expect there to be a signif- 
icant difference. We leave the investigation of clustering as 
a function of luminosity and color and exploration of results 
for higher redshifts (where the simulation bias predictions 
are much higher) to future work. 

Meneux et al (2008) have carried out a detailed study 
of clustering of galaxies as a function of stellar mass at 2 ~ 1 
(actually z = 0.85 mean redshift) in the VIMOS-VLT Deep 
Survey. This can be compared directly to our simulation 
results. To do this, we have split up our z = 1 sample 
into the same stellar mass bins as Meneux et al. . These 
are log(M/M0) = 9.5 - 10.0, log(M/M0) = 10.0 - 10.5, 
and log(M/M0) = 10.5 - 11.0. Meneux et al. find a bias 
of 1.29 ± 0.10, 1.32 ± 0.10 and 1.62 ± 0.18 for these three 
bins, respectively. Our results are 0.86 ± 0.05, 1.05 ± 0.03 
and 1.30 ± 0.07. Compared to this sample of galaxies, our 
clustering results are therefore somewhat lower. This could 
be understood if the stellar masses of our simulated galaxies 
are higher for a given space density than in the observa- 
tions. Another possibility is that the relatively small box 
size of the simulation (with the fundamental mode of the 
box going mildly non-linear at 2 = 1) is affecting the result. 
As we mentioned above, this should be less of a problem for 
the prediction of the relative bias between dark matter and 
galaxies than for the raw galaxy ^, but should still be borne 
in mind. We leave further investigation of the mass function 
by stellar mass of the simulated galaxies to future work. 

The clustering of DEEP2 galaxies around SDSS quasars 
has been measured by Coil et al. (2007), who find from the 
cross correlation that they are distributed similarly, with 



no dependence on scale, quasar luminosity or redshift. In 
the future, the BHcosmo simulation can be used to compare 
to these observations in detail. For now we note that some 
qualitatively similar findings (limited dependence of quasar 
properties on environment) has been found in the simulation 
by Colberg & Di Matteo (2008). 

6 SUMMARY AND DISCUSSION 
6.1 Summary 

We have examined the galaxies forming in a cosmological hy- 
drodynamic simulation which includes self consistent mod- 
elling of black hole growth and feedback. The mass and spa- 
tial resolution of the simulation are similar to those of previ- 
ous zoomed simulations of individual galaxies, enabling us to 
explore to a similar accuracy the properties of a population 
of ~ 1000 galaxies with circular velocities ^ 100 km s^^. 
Our conclusions can be summarised as follows. 

(1) At redshift z = 1 we find that galaxies can be roughly 
divided into classes on the basis of their morphology, kine- 
matics and age. These can be associated with the classic 
division into "early" and "late" types, so that on average 
most of the usually assumed relationships between galaxy 
properties hold. These relationships are that smaller galax- 
ies tend to be younger, have smaller Sersic indices, higher 
kinematic disk to bulge ratios and star formation rates than 
older, larger galaxies. 

(2) The angular momentum content of galaxies compared to 
observations is still a problem, as it is for the zoomed sim- 
ulations. We find that disk galaxies in the simulation have 
on average ~ 8 times less specific angular momentum for 
a given circular velocity than in the observational compi- 
lation of Navarro (1998). The kinematic decomposition of 
galaxies into disk and bulge fraction also leads to a much 
higher proportion of bulge dominated galaxies than in obser- 
vations. The inclusion of black hole physics does not appear 
to alleviate the angular problem (also found in the zoomed 
simulations of Okamoto et al. 2007). 

(3) The scale length of disk galaxies for a given circular 
velocity is approximately consistent with observations. The 
flattened parts of galaxies are therefore partially supported 
by non-circular motions. This is related to the fact that 
the simulated disk galaxies also compare acceptably to the 
z = 1 stellar mass TuUy-Fisher relation observations of 
Kassin et al. (2007), who include a non-circular component 
(gas velocity dispersion) in their preferred statistic. When 
computing the same statistic for the simulations, the overall 
agreement is acceptable, but with a different ratio of ro- 
tational and non-rotational motions than the observations. 
The extended stellar populations of galaxies are approxi- 
mately the right size, but are kinematically distributed dif- 
ferently. 

(4) The properties of galaxies in the simulation are related 
to their large scale environment, in a fashion which changes 
with redshift. At 2 = 0, the higher fraction of early types 
in dense environments is qualitatively similar to that seen 
in observational data (e.g.. Dressier 1980). At higher red- 
shifts, the observational and simulation trends are both for 
relatively more late types at high densities. However, the 
simulation predicts somewhat too many, with an actual re- 
versal of the morphology density relation. 
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(5) The overall fraction of disky (late type) galaxies in the 
simulation decreases as we move to lower redshifts, changing 
from > 0.7 at z ~ 10 to - 0.45 at 2; = 0. 

(6) The clustering of galaxies (with Vcirc > 100 km s~^) is 
slightly stronger than that of the dark matter, with the cor- 
relation function ^(r) exhibiting an average positive bias fac- 
tor b ~ 1.2 on scales r > 0.5 /i~^Mpc. At z = 1, we find 
that late type galaxies (selected by either age, kinematics or 
Sersic index) are actually more positively biased than early 
types, a reversal of the trend observed at low redshifts. 

6.2 Discussion 

As with previous zoomed simulations of disk galaxies, our 
results are mixed. There are many aspects of the simulated 
galaxies which agree qualitatively and quantitatively with 
observational data, but there are also still important differ- 
ences. The large fraction of mass in bulge components, linked 
to the lack of orbital stellar angular momentum are the most 
obvious. In order to form more realistic disk galaxies in this 
type of simulation, it is likely that subtleties in the feedback 
mechanisms need to be changed. Okamoto et al. (2005) have 
investigated different models of feedback in zoomed simula- 
tions, finding that galaxy morphologies are extremely sensi- 
tive to how the feedback is applied and how it is related to 
star formation. For example, they find that when star for- 
mation is formed through shock induced compression of gas 
(a "burst" mode), the feedback heats up a large reservoir 
of gas which later cools to form an extended rotationally 
supported disk. On the other hand, Okamoto et al. (2005) 
find that concentrating injection of feedback energy in high 
density regions (a rough approximation of AGN feedback) 
does not have the same effect. We have seen in this paper (as 
Okamoto et al. 2007 also did with zoomed simulations with 
black holes) that modelling AGN feedback with black hole 
particles does indeed not solve the angular momentum prob- 
lem. Clearly, different models of feedback should be explored 
in simulations which encompass a fair sample of galaxies. In 
the present work, we have shown that the multiphase model 
for star formation of Springel & Hernquist (2003) can be 
used to simulate galaxies in a uniform cosmological run with 
the same relative degree of success as in the zoomed runs of 
Robertson et al. (2004). 

Numerical and algorithmic issues can also infiuence the 
morphology of galaxies as well as implementations of feed- 
back physics. For example, Governato et al. (2007) have 
found that increasing numerical resolution has led to some 
improvement in the rotational support of disks and less 
transfer of angular momentum to the dark matter. We have 
seen in the present paper that our low resolution D4 run 
is not very much different that the BHCosmo run, so that 
it seems unlikely that straightforward decreases in particle 
mass and gravitational softening length would result in sig- 
nificantly bettor agreement with observations than could bo 
obtained without also changing the way star formation and 
feedback are modelled. In the present work, we have made 
use of multiphase modelling of subresolution star formation 
physics which Springel & Hernquist (2003b) have shown ex- 
hibits good convergence properties. However, better mod- 
elling of star formation and less dependence on subresolu- 
tion parametrization would go hand in hand with increased 
resolution. For example, Pelupessy et al. (2006) find that 



explicit modelling of molecular hydrogen formation in SPH 
simulations of dwarf galaxies requires gravitational soften- 
ings on the order of tens of pc and gas particle masses 
of ~ lO^h"^M0. Ceverino & Klypin (2007) have recently 
shown with AMR simulations of hot gas bubbles caused by 
star formation that modelling these bubbles self consistently 
without subresolution assumptions requires 50pc resolution. 

On the algorithmic side, there is the possiblity that is- 
sues with SPH are affecting the gas physics that can be mod- 
elled well. For example, Agertz et al. (2007) have recently 
shown that spurious pressure forces on particles can arise 
in regions where there are steep density gradients (but see 
Price, 2007 for natural ways to resolve this within the SPH 
formalism). On the other hand, Bryan (2007) has shown that 
the overly centrally concentrated density distributions and 
rotation curves that are a generic feature of simulations of 
spiral galaxies also occur in extremely high resolution sim- 
ulations of disks using Adaptive Mesh Refinement. Bryan 
(2007) concludes that more sophisticated modelling of the 
interstellar medium, including star formation and turbulent 
support is needed to resolve these issues in galaxy formation. 

Some of the results in this paper may have a slightly 
different interpretation if the classification of early and late 
type galaxies is broadened into more categories. For exam- 
ple, dwarf ellipticals and their subtype dwarf spheroidals 
are known to have many of the properties of spiral galaxies 
rather than larger ellipticals (e.g., Grebel et al. 2003). It is 
possible that some of the effects we are seeing are partly due 
to this (for example the higher clustering amplitude of late 
types at high z could partially be because there are many 
dwarf galaxies in the early type sample). Our lower cutoff 
in mass corresponds to a circular velocity ~ 100 km s^^, 
probably too high for most true dwarfs, but still low enough 
that these possiblities should be investigated further in fu- 
ture work. It is also possible that irregular galaxies and their 
properties could be investigated separately. Future work at 
higher resolution will eventually make this possible. 

The presence of late-type galaxies at high redshifts in 
high density environments which we have seen in §4.3 in 
the BHCosmo simulation has long been seen observation- 
ally. For example, clusters aX z > 0.4 contain substantial 
populations of blue spiral galax;ies (e.g., Butcher and Oem- 
ler 1978) which evolve much more rapidly than field galaxies 
in the 4-5 Gyr timescale to 2 = 0. Transformation of these 
galaxies by "galaxy harassment" (Moore et al. 1996) into 
dwarf ellipticals may be ongoing in our simulation, and par- 
tially account for the rising fraction of early types at lower 
redshifts. We can see from the morphology-density relation 
in Figure 13 that there is indeed much more evolution of the 
galaxy population at the high density end. 

Van den Bosch et al. (2007) have investigated the mor- 
phological tranformations that take place in high density en- 
vironments in the SDSS, looking at the difference between 
central and satellite galaxies. They find that whenever a cen- 
tral galax;y falls into a bigger halo (i.e., becomes a satellite), 
it undergoes a color change (becoming slightly redder) and 
a morphology change (becoming very slightly more concen- 
trated) . They find that the magnitude of this effect is inde- 
pendent of the mass of the new host halo in which the galaxy 
falls. This suggests that ram-pressure stripping of cold gas is 
not the likely cause, but instead "strangulation" , or removal 
of the galaxy halo hot gas. This takes place over a relatively 
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long timescale (not accounted for in past semi-analytic mod- 
els, as shown by Kang & van den Bosch (2008). In our sim- 
ulation, we do find in Figure 12 that galaxies become more 
centrally concentrated (n increases) function of de- 

creasing distance from the host halo center. Tracking the 
role of gas accretion and loss from the galaxies (as carried 
out by e.g., Keres et al. 2005) could be used to investigate 
the mechanism by which this occurs. 

Much of the focus of galaxy formation simulations has 
been on trying to form realistic disks. We have seen in the 
present work that we are able to achieve a similar degree 
of success (for example in the specific angular momentum- 
circular velocity relationship) in our full volume cosmologi- 
cal run to earlier zoomed simulations (e.g., Robertson et al. 
2004, Abadi et al. 2003ab). We have not looked in so much 
detail however at the elliptical galaxies formed in the BH- 
Cosmo. Naab et al. (2007) have recently performed zoomed 
simulations of early-type galaxies from ACDM initial and 
find that realistic intermediate mass elliptical galaxies can 
be formed even without recent major mergers or AGN and 
supernova feedback. We have mass resolution throughout 
our box only ~ 3 times worse, than in the fiducial zoomed 
runs of Naab et al. (although we have significantly coarser 
spatial resolution), so that it would be worth investigating 
the properties of the larger population of early types in more 
detail in future work. 

In the future, it will be necessary to run a larger box 
down to redshift z = 0. \n this way, it will be possible to 
investigate a wider range of environments, including rich 
clusters, where the evolutionary history of simulated galax- 
ies is expected to be different from galaxy groups and the 
field. A much wider range of observational data is also avail- 
able at 2; = 0. Population synthesis modelling of the sim- 
ulated galaxies would enable more direct comparisons with 
the data, including the effect of observing in different wave- 
lengths on morphology. 
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